﻿// OCLSmallMatrices.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "OpenCLKernel.h"
#include <math.h>
#include <time.h>


const char* _kernel_multiply_00_ = KERNEL(
    __kernel void multiply_00(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
    __local float4 Asub[4], Bsub[4];
    const int li0 = get_local_id(0), li1 = get_local_id(1);
    const int gr0 = get_group_id(0), gr1 = get_group_id(1);
    float acc = 0.0f;

    Asub[li0] = vload4(((gr0 * 4 + li0) * colA) / 4, inputA);
    Bsub[li1] = (float4)(inputB[gr1 * 4 + li1], inputB[colB + gr1 * 4 + li1], inputB[2 * colB + gr1 * 4 + li1], inputB[3 * colB + gr1 * 4 + li1]);
    //barrier(CLK_LOCAL_MEM_FENCE);

    acc = dot(Asub[li0], Bsub[li1]);
    //barrier(CLK_LOCAL_MEM_FENCE);

    outputC[(gr0 * 4 + li0) * colB + gr1 * 4 + li1] = acc;

    return;
}
);


const char* _kernel_multiply_01_ = KERNEL(
    __kernel void multiply_01(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
    const int i = get_global_id(0), j = get_global_id(1);
    int k;
    float sum = 0.0f;
    for (k = 0; k < colA; k++)
        sum += inputA[i * colA + k] * inputB[k * colB + j];
    outputC[i * colB + j] = sum;
    return;
}
);


// merge 所有小矩阵，一起运算（注意16即为Kd, inputA, inputB, inputC的尺寸为所有小矩阵）
const char* _kernel_multiply_02_ = KERNEL(
    __kernel void multiply_02(__global float* inputA, __global float* inputB, __global float* outputC, int rowA, int colA, int colB)
{
    __local float4 Asub[4][16], Bsub[4][16];
    const int li0 = get_local_id(0), li1 = get_local_id(1);
    const int gr0 = get_group_id(0), gr1 = get_group_id(1);
    float acc[16] = { 0.0f };
    int i;

    for (i = 0; i < 16; i++)
    {
        Asub[li0][i] = vload4(((gr0 * 4 + li0) * colA * 16 + i * 4) / 4, inputA);
        Bsub[li1][i] = vload4(((gr1 * 4 + li1) * colA * 16 + i * 4) / 4, inputB);
    }
    //barrier(CLK_LOCAL_MEM_FENCE);

    for (i = 0; i < 16; i++)
        acc[i] = dot(Asub[li0][i], Bsub[li1][i]);
    //barrier(CLK_LOCAL_MEM_FENCE);

    for (i = 0; i < 16; i++)
        outputC[i * rowA * colB + (gr0 * 4 + li0) * colB + gr1 * 4 + li1] = acc[i];
    return;
}
);


const int rowA = 1600, colA = 4, colB = 800;
const int Kd = 4 * 4; // K 路 d 阶个小矩阵相乘


bool floatEq(const float a, const float b)// 相等返回 1
{
    if (b == 0)
        return fabs(a) < 0.001;
    return fabs(1.0 * a / b - 1) < 0.001;
}


// 顺序执行Kd个
long SerialProcess(OpenCLKernel* pOCL, cl_kernel* l_kernel, cl_mem* l_bufferA, cl_mem* l_bufferB, cl_mem* l_bufferC, float* outC)
{
    cl_int status;

    cl_command_queue queue = clCreateCommandQueue(pOCL->GetContext(), pOCL->GetDevices()[0], NULL, &status);
    GET_CL_ERROR(status, pOCL);

    clock_t time_start = clock();

    //参数赋值
    for (int i = 0; i < Kd; i++)
    {
        status = clSetKernelArg(l_kernel[i], 0, sizeof(cl_mem), &l_bufferA[i]);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 1, sizeof(cl_mem), &l_bufferB[i]);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 2, sizeof(cl_mem), &l_bufferC[i]);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 3, sizeof(int), &colA);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 4, sizeof(int), &colB);
        GET_CL_ERROR(status, pOCL);
    }

    //工作组
    size_t* globalSize = NULL;
    size_t* localSize = NULL;
    cl_uint work_dim = 0;
    work_dim = 2;
    globalSize = new size_t[2];
    localSize = new size_t[2];
    globalSize[0] = rowA;
    globalSize[1] = colB;
    localSize[0] = 4;
    localSize[1] = 4;

    //执行
    for (int i = 0; i < Kd; i++)
    {
        status = clEnqueueNDRangeKernel(queue, l_kernel[i], work_dim, NULL, globalSize, localSize, 0, NULL, NULL);
        GET_CL_ERROR(status, pOCL);
        printf("%d\n", i);
    }
    status = clFinish(queue);
    GET_CL_ERROR(status, pOCL);

    clock_t time_end = clock();

    //写内存
    status = clEnqueueReadBuffer(queue, l_bufferC[0], CL_TRUE, 0, sizeof(float) * rowA * colB, outC, 0, NULL, NULL);
    GET_CL_ERROR(status, pOCL);

    if (globalSize != NULL)
        delete[] globalSize;
    if (localSize != NULL)
        delete[] localSize;
    clReleaseCommandQueue(queue);

    return time_end - time_start;
}


// 异步执行Kd个
long AsynchProcess(OpenCLKernel* pOCL, cl_kernel* l_kernel, cl_mem* l_bufferA, cl_mem* l_bufferB, cl_mem* l_bufferC, float* outC)
{
    cl_int status;

    cl_command_queue l_queue[Kd];
    for (int i = 0; i < Kd; i++)
    {
        l_queue[i] = clCreateCommandQueue(pOCL->GetContext(), pOCL->GetDevices()[0], NULL, &status);
        GET_CL_ERROR(status, pOCL);
    }

    cl_event GPUExecution[Kd];

    clock_t time_start = clock();

    //参数赋值
    for (int i = 0; i < Kd; i++)
    {
        status = clSetKernelArg(l_kernel[i], 0, sizeof(cl_mem), &l_bufferA[i]);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 1, sizeof(cl_mem), &l_bufferB[i]);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 2, sizeof(cl_mem), &l_bufferC[i]);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 3, sizeof(int), &colA);
        GET_CL_ERROR(status, pOCL);
        status = clSetKernelArg(l_kernel[i], 4, sizeof(int), &colB);
        GET_CL_ERROR(status, pOCL);
    }

    //工作组
    size_t* globalSize = NULL;
    size_t* localSize = NULL;
    cl_uint work_dim = 0;
    work_dim = 2;
    globalSize = new size_t[2];
    localSize = new size_t[2];
    globalSize[0] = rowA;
    globalSize[1] = colB;
    localSize[0] = 4;
    localSize[1] = 4;

    //执行
    for (int i = 0; i < Kd; i++)
    {
        status = clEnqueueNDRangeKernel(l_queue[i], l_kernel[i], work_dim, NULL, globalSize, localSize, 0, NULL, &GPUExecution[i]);
        GET_CL_ERROR(status, pOCL);
        printf("%d\n", i);
    }
    status = clWaitForEvents(Kd, GPUExecution);
    GET_CL_ERROR(status, pOCL);

    clock_t time_end = clock();

    //写内存
    status = clEnqueueReadBuffer(l_queue[0], l_bufferC[0], CL_TRUE, 0, sizeof(float) * rowA * colB, outC, 0, NULL, NULL);
    GET_CL_ERROR(status, pOCL);

    if (globalSize != NULL)
        delete[] globalSize;
    if (localSize != NULL)
        delete[] localSize;
    for (int i = 0; i < Kd; i++)
        clReleaseCommandQueue(l_queue[i]);

    return time_end - time_start;
}


// 合并执行
long MergeProcess(OpenCLKernel* pOCL, cl_kernel kernel, cl_mem bufferA, cl_mem bufferB, cl_mem bufferC, float* outC)
{
    cl_int status;

    cl_command_queue queue = clCreateCommandQueue(pOCL->GetContext(), pOCL->GetDevices()[0], NULL, &status);
    GET_CL_ERROR(status, pOCL);

    clock_t time_start = clock();

    //参数赋值
    status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufferA);
    GET_CL_ERROR(status, pOCL);
    status = clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufferB);
    GET_CL_ERROR(status, pOCL);
    status = clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufferC);
    GET_CL_ERROR(status, pOCL);
    status = clSetKernelArg(kernel, 3, sizeof(int), &rowA);
    GET_CL_ERROR(status, pOCL);
    status = clSetKernelArg(kernel, 4, sizeof(int), &colA);
    GET_CL_ERROR(status, pOCL);
    status = clSetKernelArg(kernel, 5, sizeof(int), &colB);
    GET_CL_ERROR(status, pOCL);

    //工作组
    size_t* globalSize = NULL;
    size_t* localSize = NULL;
    cl_uint work_dim = 0;
    work_dim = 2;
    globalSize = new size_t[2];
    localSize = new size_t[2];
    globalSize[0] = rowA;
    globalSize[1] = colB;
    localSize[0] = 4;
    localSize[1] = 4;

    //执行
    status = clEnqueueNDRangeKernel(queue, kernel, work_dim, NULL, globalSize, localSize, 0, NULL, NULL);
    GET_CL_ERROR(status, pOCL);
    status = clFinish(queue);
    GET_CL_ERROR(status, pOCL);

    clock_t time_end = clock();

    //写内存
    status = clEnqueueReadBuffer(queue, bufferC, CL_TRUE, 0, sizeof(float) * rowA * colB, outC, 0, NULL, NULL);
    GET_CL_ERROR(status, pOCL);

    if (globalSize != NULL)
        delete[] globalSize;
    if (localSize != NULL)
        delete[] localSize;
    clReleaseCommandQueue(queue);

    return time_end - time_start;
}


int main()
{
    float* A = (float*)malloc(sizeof(float) * rowA * colA);
    float* B = (float*)malloc(sizeof(float) * colA * colB);
    float* C = (float*)malloc(sizeof(float) * rowA * colB);
    srand(2);
    for (int i = 0; i < rowA * colA; A[i] = rand() & 0xF, i++);
    for (int i = 0; i < colA * colB; B[i] = rand() & 0xF, i++);
    for (int i = 0; i < rowA * colB; C[i] = 0, i++);

    cl_int status;

    OpenCLKernel* m_pOCL = new OpenCLKernel;
    if (0 != m_pOCL->OpenCLInitialize())
    {
        delete m_pOCL;
        m_pOCL = nullptr;
        return EXIT_FAILURE;
    }

    //cl_command_queue queue = clCreateCommandQueue(m_pOCL->GetContext(), m_pOCL->GetDevices()[0], NULL, &status);
    //GET_CL_ERROR(status, m_pOCL);

    //编译
    size_t nSourceSize = NULL;
    cl_program program = NULL;
    //cl_kernel kernel = NULL;

    // 0.多kernel顺序;1.多kernel异步;2.单kernel批量
    printf("Choose the way of multiply.\n");
    int sel = 0;
    scanf_s("%d", &sel, 1);

    if (sel == 0 || sel == 1)
    {
        cl_kernel l_kernel[Kd];

        nSourceSize = { strlen(_kernel_multiply_00_) };
        program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_multiply_00_, &nSourceSize, &status);
        GET_CL_ERROR(status, m_pOCL);
        status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
        GET_CL_ERROR(status, m_pOCL);
        for (int i = 0; i < Kd; i++)
        {
            l_kernel[i] = clCreateKernel(program, "multiply_00", &status);
            GET_CL_ERROR(status, m_pOCL);
        }

        //读内存
        cl_mem l_bufferA[Kd], l_bufferB[Kd], l_bufferC[Kd];
        for (int i = 0; i < Kd; i++)
        {
            l_bufferA[i] = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * rowA * colA, A, &status);
            GET_CL_ERROR(status, m_pOCL);
            l_bufferB[i] = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * colA * colB, B, &status);
            GET_CL_ERROR(status, m_pOCL);
            l_bufferC[i] = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_WRITE_ONLY, sizeof(float) * rowA * colB, NULL, &status);
            GET_CL_ERROR(status, m_pOCL);
        }

        clock_t proc_time = 0;
        if (sel == 0)
            proc_time = SerialProcess(m_pOCL, l_kernel, l_bufferA, l_bufferB, l_bufferC, C);
        else
            proc_time = AsynchProcess(m_pOCL, l_kernel, l_bufferA, l_bufferB, l_bufferC, C);
        printf("\nTime kernel : %d ms\n", proc_time);

        for (int i = 0; i < Kd; i++)
        {
            clReleaseMemObject(l_bufferA[i]);
            clReleaseMemObject(l_bufferB[i]);
            clReleaseMemObject(l_bufferC[i]);
        }
        for (int i = 0; i < Kd; i++)
            clReleaseKernel(l_kernel[i]);
    }
    else if (sel == 2)
    {
        cl_kernel kernel = NULL;

        nSourceSize = { strlen(_kernel_multiply_02_) };
        program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_multiply_02_, &nSourceSize, &status);
        GET_CL_ERROR(status, m_pOCL);
        status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
        GET_CL_ERROR(status, m_pOCL);
        kernel = clCreateKernel(program, "multiply_02", &status);
        GET_CL_ERROR(status, m_pOCL);

        //读内存
        cl_mem bufferA, bufferB, bufferC;

        float* new_A = (float*)malloc(sizeof(float) * rowA * colA * Kd);
        for (int j = 0; j < rowA; j++)
            for (int i = 0; i < Kd; i++)
                memcpy(new_A + j * colA * Kd + i * colA, A + j * colA, sizeof(float) * colA);
        bufferA = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * rowA * colA * Kd, new_A, &status);
        GET_CL_ERROR(status, m_pOCL);

        float* B_T = (float*)malloc(sizeof(float) * colA * colB);
        for (int x = 0; x < colA; x++)
            for (int y = 0; y < colB; y++)
                B_T[y * colA + x] = B[x * colB + y];
        float* new_B = (float*)malloc(sizeof(float) * colA * colB * Kd);
        for (int j = 0; j < colB; j++)
            for (int i = 0; i < Kd; i++)
                memcpy(new_B + j * colA * Kd + i * colA, B_T + j * colA, sizeof(float) * colA);
        bufferB = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * colA * colB * Kd, new_B, &status);
        GET_CL_ERROR(status, m_pOCL);

        bufferC = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_WRITE_ONLY, sizeof(float) * rowA * colB * Kd, NULL, &status);
        GET_CL_ERROR(status, m_pOCL);

        clock_t proc_time = 0;
        proc_time = MergeProcess(m_pOCL, kernel, bufferA, bufferB, bufferC, C);
        printf("\nTime kernel : %d ms\n", proc_time);

        free(new_A);
        free(B_T);
        free(new_B);
        clReleaseMemObject(bufferA);
        clReleaseMemObject(bufferB);
        clReleaseMemObject(bufferC);
        clReleaseKernel(kernel);
    }

    // 返回并检查结果
    printf("Checking...\n");
    int correct = 1;
    for (int i = 0; i < rowA && correct; i++)
    {
        for (int j = 0; j < colB && correct; j++)
        {
            float tempSum = 0.0f;
            for (int k = 0; k < colA; k++)
            {
                tempSum += A[i * colA + k] * B[k * colB + j];
            }
            if (!floatEq(tempSum, C[i * colB + j]))
            {
                printf("Error at [%d, %d], calculation: %f, reference: %f\n", i, j, C[i * colB + j], tempSum);
                correct = 0;
            }
        }
    }
    if (correct)
        printf("Result correct.\n");

    //结束并释放
    clReleaseProgram(program);
    delete m_pOCL;
    free(A);
    free(B);
    free(C);

    return 0;
}
